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ABSTRACT 

We develop a new symbolic-numeric algorithm for the certi- 
fication of singular isolated points, using their associated lo- 
cal ring structure and certified numerical computations. An 
improvement of an existing method to compute inverse sys- 
tems is presented, which avoids redundant computation and 
reduces the size of the intermediate linear systems to solve. 
We derive a one-step deflation technique, from the descrip- 
tion of the multiplicity structure in terms of differentials. 
The deflated system can be used in Newton-based iterative 
schemes with quadratic convergence. Starting from a poly- 
nomial system and a small-enough neighborhood, we obtain 
a criterion for the existence and uniqueness of a singular root 
of a given multiplicity structure, applying a well-chosen sym- 
bolic perturbation. Standard verification methods, based 
eg. on interval arithmetic and a fixed point theorem, are 
employed to certify that there exists a unique perturbed 
system with a singular root in the domain. Applications to 
topological degree computation and to the analysis of real 
branches of an implicit curve illustrate the method. 

Categories and Subject Descriptors 

G.1.5 [Mathematics of Computing]: Roots of Nonlinear 
Equations; 1.1.2 [Computing Methodologies]: Symbolic 
and Algebraic Manipulation — Algebraic algorithms 

General Terms 

Algorithms, Theory 

Keywords 

root deflation, multiplicity structure, dual space, inverse sys- 
tem, isolated point 

1. INTRODUCTION 

A main challenge in algebraic and geometric computing 
is singular point identification and treatment. Such prob- 
lems naturally occur when computing the topology of im- 
plicit curves or surfaces pTJ, the intersection of parametric 
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surfaces in geometric modeling. When algebraic representa- 
tions are used, this reduces to solving polynomial systems. 
Several approaches are available: algebraic techniques such 
as Grobner bases or border bases, resultants, subdivision al- 
gorithms[10|. |13j . homotopies, and so on. At the end of 
the day, a numerical approximation or a box of isolation is 
usually computed to characterize the (real) roots of the poly- 
nomial system. But we often need to improve the numerical 
approximation of the roots. Numerical methods such that 
Newton's iteration can be used to improve the quality of 
the approximation, provided that we have a simple root. In 
the presence of a multiple root, the difficulties are signifi- 
cantly increasing. The numerical approximation can be of 
very bad quality, and the methods used to compute this ap- 
proximation are converging slowly (or not converging). The 
situation in practical problems, as encountered in CAGD 
for instance, is even worse, since the coefficients of the in- 
put equations are known, with some incertitude. Computing 
multiple roots of approximate polynomial systems seems to 
be an ill-posed problem, since changing slightly the coeffi- 
cients may transform a multiple root into a cluster of simple 
roots (or even make it disappear). 

To tackle this difficult problem, we adopt the following 
strategy. We try to find a (small) perturbation of the input 
system such that the root we compute is an exact multi- 
ple root of this perturbed system. In this way, we identify 
the multiplicity structure and we are able to setup deflation 
techniques which restore the quadratic convergence of the 
Newton system. The certification of the multiple root is also 
possible on the symbolically perturbed system by applying 
a fixed point theorem, based eg. on interval arithmetic [16] 
or a-theorems |17| . 

Related work. In order to be develop Newton-type meth- 
ods that converge to multiple roots, deflation techniques 
which consist in adding new equations in order to reduce 
the multiplicity have already been considered. In [14], by 
applying a triangulation preprocessing step on the Jacobian 
matrix at the approximate root, minors of the Jacobian ma- 
trix are added to the system to reduce the multiplicity. 

In [BJ, a presentation of the ideal in a triangular form in 
a good position and derivations with respect to the leading 
variables are used to iteratively reduce the multiplicity. This 
process is applied for p-adic lifting with exact computation. 

In [8], instead of triangulating the Jacobian matrix, 
the number of variables is doubled and new equations are 
introduced, which are linear in the new variables. They des- 
cribe the kernel of the Jacobian matrix at the multiple root. 

In [3], this construction in related to the construction of 



the inverse system. The dialytic method of F.S. Macaulay 
[9] is revisited for this purpose. These deflation methods are 
applied iteratively until the root becomes regular, doubling 
each time the number of variables. 

More recent algorithms for the construction of inverse sys- 
tems are described eg. in [11], reducing the size of the in- 
termediate linear systems (and exploited in [TS]), ° r in [12] 
using an integration method. 

In [15], a minimization approach is used to reduce the 
value of the equations and their derivatives at the approxi- 
mate root, assuming a basis of the inverse system is known. 

In [2D], the inverse system is constructed via Macaulay 's 
method; tables of multiplications are deduced and their eigen- 
values are used to improve the approximated root. It is 
proved that the convergence is quadratic when the Jacobian 
has corank 1 at the multiple root. 

Verification of multiple roots of (approximate) polynomial 
equations is a difficult task. The approach proposed in [16] 
consists in introducing perturbation parameters and to cer- 
tifying the multiple root of nearby system by using a fixed 
point theorem, based on interval arithmetic. It applies only 
to cases where the Jacobian has corank equal to 1. 
The univariate case. In preparation for the multivariate 
case, we review some techniques used to treat singular zeros 
of univariate polynomials, and we present our method on a 
univariate instance. 

Let g(x) G K[x] be a polynomial which attains at x — 
a root of multiplicity n > 1. The latter is defined as the 
smallest positive integer fi such that d^g(0) 7^ whereas 
3(0) = dg(0) = ••• = d M_1 p(0) = 0. Here we denote 
d k g(x) = -j-jzgix) /k\ the normalized k— th order derivative. 

We see that &o = (1, d 1 , . . . , is the maximal space 

of differentials which is stable under derivation, that vanish 
when applied to members of Qo, the (x)— primary compo- 
nent of (g) at x — 0. Consider now the symbolically per- 
turbed equation 



fi(x, e) = g(x) + £1 + e 2 x H h£ M -2£ M (1) 

and apply every basis element of Do to arrive to the new sys- 
tem f(x,e)— (fi,difi,...,d^~ 1 fi \ in n — 1 variables. The 

i-th equation shall be = cT^fi = d^g+^Zi x k ~ i+1 e k , 
i.e linear in e, the last one being / M = d^~ g(x). This sys- 
tem deflates the root, as we see that the determinant of it's 
Jacobian matrix at (0, 0) is 



det J f (0, 0) = 



d f 
d^J 1 


1 


d f 


1 








-Mm (0,0) 
-IMd"g{0) ^ 0. 



Now suppose that f* is an approximate zero, close to x = f . 
We can still compute 2\ by evaluating g(x) and the deriva- 
tives up to a threshold relative to the error in £*. Then 
we can form ([l| and use verification techniques to certify 
the root. Checking that the Newton operator is contract- 
ing shows the existence and unicity of a multiple root in a 
neighborhood of the input data. We are going to extend 
this approach, described in [15], to multi-dimensional iso- 
lated multiple roots. 

Our approach. It consists of the following steps: 

(a) Compute a basis for the dual space and of the local 
quotient ring at a given (approximate) singular point. 



(b) Deflate the system by augmenting it with new equa- 
tions derived from the dual basis, introducing adequate per- 
turbation terms. 

(c) Certify the singular point and its multiplicity structure 
for the perturbed system checking the contraction property 
of Newton iteration (eg. via interval arithmetic). 

In step (a), a dual basis at the singular point is computed 
by means of linear algebra, based on the integration ap- 
proach of [15]. We describe an improvement of this method, 
which yields directly a triangular dual basis with no redun- 
dant computation. This method has the advantage to re- 
duce significantly the size of the linear systems to solve at 
each step, compared to Macaulay's type methods [5] [7] [S] 
[3]. In the case of an approximate singular point, errors are 
introduced in the coefficients of the basis elements. Yet a 
successful computation is feasible. In particular, the sup- 
port of the basis elements is revealed by this approximate 
process. 

In the deflation step (b), new equations and new variables 
are introduced in order to arrive to a new polynomial system 
where the singularity is obviated. The new variables corre- 
spond to perturbations of the initial equations along specific 
polynomials, which form a dual counterpart to the basis of 
the dual space. One of the deflated systems that we compute 
from the dual system is a square n x n system with a sim- 
ple root. This improves the deflation techniques described 
in |?1 |8l |3], which require additional variables and possi- 
bly several deflation steps. New variables are introduced 
only in the case where we want to certify the multiplicity 
structure. The perturbation techniques that we use extend 
the approach of [16] to general cases where the co-rank of 
the Jacobian matrix could be bigger than one. The verifi- 
cation step (c) is mostly a contraction condition, using eg. 
techniques as in [16]. This step acts on the (approximate) 
deflated system, since verifying a simple solution of the de- 
flated system induces a certificate of an exact singular point 
of (a nearby to) the initial system. 

We are going to detail the different steps in the following 
sections, starting with notations in Sect. [5] dual basis in 
Sect. [31 deflation in Sect. [3] and certification in Sect. [S] In 
the last section, we will show examples and applications to 
the topology analysis of curves. 

2. PRELIMINARIES AND MAIN RESULTS 

We denote by R = K[x] a polynomial ring over the field 
K of characteristic zero. Also, the dual ring R* is the space 
of linear functionals A : R — ¥ K. It is commonly identified 
as the space of formal series K[[e?]] where d = (<9i, . . . , <9„) 
are formal variables. Thus we view dual elements as formal 
series in differential operators at a point £ G K" . To specify 
that we use the point £, we also denote these differentials d^. 
When applying A(d() G K[[9^]] to a polynomial g(x) G R 
we will denote by A**[g] = A^g — A(d^)[g(x)] the operation 



Oil'. 



d? 1 • • • d? 



-(0, 



(2) 



for A(d c ) = G K[[0 C ]]. Extending this definition 

to an ordered set T> = (Ai, . . . , A M ) G K[[9]] M , we shall de- 
note T)^ [g] = (A^gr, . . . , A^g). In some cases, it is convenient 
to use normalized differentials instead of d: for any a G N n , 
we denote d% = i d%. When C, = 0, we have d^x* 3 = 1 



if a — (3 and otherwise. More generally, (dJ?) aS pjn is the 
dual basis of ((a; — £) a ) ae n™- 

For A G R* and p G R, let p ■ A : q M> A(pq). We check 
that 

( ari _C0-o? = ^(af). (3) 

This property shall be useful in the sequel. 

2.1 Isolated points and differentials 

Let (/) = {f%, ... , f s ) be an ideal of R, £ G K a root of 
/ and m<; = (xi — fi, . . . ,x n — C«) the maximal ideal at £. 
Suppose that £ is an isolated root of /, then a minimal pri- 
mary decomposition of X = O Q contains a primary 

S prim. 31 



component Q^ such that y Q<; = and Vy m( for the 
other primary components Q' associated to X [2]. 

As i/Oc = m C> R/Q-C IS a finite dimensional vector space. 
The multiplicity p^ of £ is defined as the dimension of R/Qc- 
A point of multiplicity one is called regular point, or simple 
root, otherwise we say that £ is a singular isolated point, or 
multiple root of /. In the latter case we have Jf (C) = 0. 

We can now define the dual space of an isolated point. 

Definition 2.1. The dual space of X is the subspace of 
elements ofK[[d(]] that vanish on all the elements ofX. It 
is also called the orthogonal ofX and denoted by X ± . 

Consider now the orthogonal of Qc, i.e. the subspace &c 
of elements of R* that vanish on members of Qc, namely 

Qc = & c = {A G R* : A c [p] = 0, Vp G Q c }. 

The following lemma is an essential property that allows ex- 
traction of the local structure S>c directly from the "global" 
ideal X = (/) , notably by matrix methods outlined in Sect. [3] 



Proposition 2.2 (Q21 Th. 8]). For any isolated point 
CeKof f, we have (/) ± n K[d c ] = S> c . 

In other words, we can identify &c = with the space of 
polynomial differential operators that vanish at £ on every 
element of X. 

The space Q>c is a vector space of polynomials in dc di- 
mension the multiplicity of £. As the variables (xi — £j) 
act on R* as derivations (see (O), 2#c is a space of differen- 
tial polynomials in 8q, which is stable by derivation. This 
property will be used explicitly in constructing Sc. (Sec. [3]). 

Definition 2.3. The nilindex of Qc is the maximal in- 
teger N G N s.t m£ f- Qc- 

It is directly seen that the maximal degree of an element of 
2$C has to be equal to N , also known as the depth of &c- 

2.2 Quotient ring and dual structure 

In this section we explore the relation between the dual 
ring and the quotient R/ Qc, where Qc is the primary com- 
ponent of the isolated point £. We show how to extract a 
basis of this quotient ring from the support of the elements 
of £#c an d how 2>c can be used to reduce any polynomial 
modulo Qc- 

It is convenient in terms of notation to make the assump- 
tion C = 0. This saves some indices, while it poses no con- 
straint (since it implies a linear change of coordinates), and 
shall be adopted hereafter and in the next section. 



Let suppi^o be the set of exponents of monomials appear- 
ing in *2>q, with a non-zero coefficient. These are of degree 
at most N, the nilindex of Q . Since VA £ f , A°[p] = iff 
p G Qo, we derive that supp^o = {a : x a (/ Qo}- In par- 
ticular, we can find a basis of Rj Qo between the monomials 
{x a : a G supp£^}. This is a finite set of monomials, since 
their degree is bounded by the nilindex of Qo- 

Given a standard basis B = (xP i )j=i,..., M of R/Qo and, for 
all monomials x^i <f Qo, j = 1, ■ ■ ■ , s — p, s — #supp £f , 
with x' 1 ' ^ B, the expression (normal form) of 

x~<i = X i:j x^ mod Qo (4) 

i=l 

in the basis B then the dual elements |12l Prop. 13] 

A< = <*"*+ ^ AucT', (5) 

j=i 

for i = 1, . . . , p, form a basis of S>o- We give a proof of this 
fact in the following lemma. 

Lemma 2.4. The set of elements T> = (Aj)^!,...,^ is a ba- 
sis of &o and the normal form of any g(x) G R with respect 



to the standard basis B = (x l ^ i )i = 
AF(.g) = f>°[ 



(6) 



Proof. First note that the elements of T> axe linearly in- 
dependent. Now, by construction, I]f =1 A°[a; a ] = NF(a; a ) 
for all x a i Qo, eg. NF(cc /3 *) = x Pi . Also, for x a G Qo, 
Vi, A; (x a ) = 0, since a (£ suppO. Thus the elements of T> 
compute NF(-) on all monomials of R, and ([6]) follows by li- 
nearity. We deduce that T> is a basis of S)o, as in Def. 12.11 □ 

Computing the normal form of the border monomials of B 
via © also yields the border basis relations and the opera- 
tors of multiplication in the quotient R/Qo (see eg. [5] for 
more properties). 

If a graded monomial ordering is fixed and B = (xP i )i=i,.., M 
is the corresponding standard basis of R/Qo, then dP 1 is the 
leading term of ([5| wrt the opposite ordering [8] Th. 3.1]. 

Conversely, if we are given a basis T> of $)o whose co- 
efficient matrix in the dual monomials basis {d a ) a ^Q is 
D G K MXS , we can compute a basis of R/Qo by choosing p, 
independent columns of D, say those indexed by dP i , i = 
1, . . . , n . If G G K MXM is the (invertible) matrix formed by 
these columns, then D' := G~ 1 D, is 



D' 



A' 



Pi 
1 





0, 




7i 
Ai,i 



1 A 



7., 
Ai. 



(7) 



i.e. a basis of the form Note that an arbitrary basis of & 
does not have the above diagonal form, nor does it directly 
provide a basis for R/Qo- 

For 1 £ N, ^ denotes the vector space of polynomials of 
2> of degree < t. The Hilbert function h : N — > N is defined 
by h(t) = dim(^t), t > 0, hence h(0) = 1 and h(t) = dim& 
for t > N. The integer h(l) — 1 = corank Jf is known as the 
breadth of 



3. COMPUTING LOCAL RING STRUCTURE 

The computation of a local basis, given a system and a 
point, is done essentially by matrix- kernel computations, 
and consequently it can be carried out numerically, even 
when the point or even the system is inexact. Throughout 
the section we suppose / £ R m and £ £ K* 1 with /(£) = 0. 

Several matrix constructions have been proposed, that use 
different conditions to identify the dual space as a null-space. 
They are based on the stability property of the dual basis: 



VAef,, ^Aefw 

ddi 



,n. 



(8) 



We list existing algorithms that compute dual-space bases: 

• As pointed out in ([3]), an equivalent form of ([8]) is: VA £ 
@t,A[xifj] = 0, Vi, j = 1, . . . , n. Macaulay's method [5] uses 
it to derive the algorithm that is outlined in Sect. 13.11 

• In [IB they exploit © by forming the matrix Di of the 

map : K[9]t — > K[d]t-i for all i = 1, . . . , n and some 
ddi 

triangular decomposition of the differential polynomials in 
terms the differential variables. This approach was used in 
|18| to reduce the row dimension of Macaulay's matrix, but 
not the column dimension. The closedness condition is also 
used in [21] to identify a superset of suppii?t+i. 

• The integration method in [T5] "integrates" elements of 
a basis of &t, and obtains a priori knowledge of the form of 
elements in degree t+1 (Sect.[3]2}. 

All methods are incremental, in the sense that they start 
by setting T>o ~ (1) and continue by computing T>i, i = 
1, . . . , N, N + 1. When #V N = #V N+1 then V N is a basis 
of &, and N is the nilindex of Q. 

We shall review two of these approaches to compute a 
basis for @, and then describe an improvement, that al- 
lows simultaneous computation of a quotient ring basis while 
avoiding redundant computations. 

3.1 Macaulay's dialytic matrices 

This matrix construction is presented in [5] Ch. 4], a mo- 
dern introduction is contained in [3], together with an im- 
plementation of the method in ApaTool^j. 

The idea behind the algorithm is the following: An ele- 
ment of Q is of the form A(&) = \ a d a under the 

\a\<N 

condition: A evaluates to at any g £ (/}, i.e. A°(<;) = 



A°(E.9*/0 = 



A°(x l3 fi) = for all monomials x^ 



If we apply this condition recursively for |a| < N we get 
a vector of coefficients (A a )| a |<jv in the (right) kernel of 
the matrix with rows indexed by constraints A°[x (3 f i ] = 0, 
\/3\<N-l. 

Note that the only requirement is to be able to perform 
derivation of the input equations and evaluation at £ = 0. 

Example 3.1. Let fi = x\— X2+x\, ji = xi— Xi+x%. We 
also refer the reader to [3] Ex. 2] for a detailed demonstration 
of Macaulay's method on the same instance. The matrices 
in order 1 and 2 are: 



fi 
h 



d 2 
-1 
-1 





1 


di 


d> 


di 


did2 


4 


h 


"0 


1 


-1 


1 








h 





1 


-1 








1 


XI fl 











1 


-1 
















1 


-1 



















1 


-1 


Z2/2 














1 


-1 



The kernel of the left matrix gives T>\ = (1, di + ^2). Ex- 
panding up to order two, we get the matrix on the right, 
and T>2 = (1, di + di- —d\ + d'{ + d\d2 + If we expand 
up to depth 3 we get the same null-space, thus T> — T>2. □ 

3.2 Integration method 

This method is presented in [12]. It is an evolution of 
Macaulay's method, in the sense that the matrices are not 
indexed by all differentials, but just by elements based on 
knowledge of the previous step. This performs a computa- 
tion adapted to the given input and results in smaller ma- 
trices. 

For A £ K[&], we denote by J k A the element $ £ K[d] 
with the property j§-&(d) = A(t?) and with no constant 
term wrt dk- 



Theorem 3.2 ( [H Th. 15] ). Let (Ai, A 2 , . . . , A s ) be 

a basis of 2$t-i, that is, the subspace of & of elements of or- 
der at most t — 1. An element A £ ¥L[d] with no constant 
term lies in Qlt iff it is of the form: 

S 71 

A ( a ) = EEW fc A *(di,...,c\,0,...,0), (9) 

i=l k = l 

for \it £ K, and the following two conditions hold: 



d 



for all 1 < k < I < n . 

(11) A<[f k ] =0, for fc = l,...,n . 

Condition (i) is equivalent to ^g- A £ S't-i, for 1 < k < n. 
Thus the two conditions express exactly the fact that *3 
must be stable by derivation and his members must vanish 
on </). 

This gives the following algorithm to compute the dual 
basis: Start with T>o — (1). Given a basis of S>t-i we gener- 
ate the ns candidate elements J k Ai_i(ft, . . . , dk, 0, . . . , 0). 
Conditions (i) and (ii) give a linear system with unknowns 
\ik. The columns of the corresponding matrix are indexed 
by the candidate elements. Then the kernel of the matrix 
gives the possible solutions, which we add to T>t-i to obtain 
T>t. If for some t there are no further new elements, then 
T> = T>t is a basis of &. 

Example 3.3. Consider the instance of Ex. 13. II We have 
/i(C) = /2(C) = 0, thus we set V = {1}. Equation © gives 
A = Aidi + \2d2- Condition (i) induces no constraints and 
(ii) yields the system 



1 -1 




" Ai ' 


1 -1 




A 2 







(10) 



where the columns are indexed by d\,d2- We get Ai = A2 = 
1 from the kernel of this matrix, thus T>\ — {1, d\ +6(2}- 

For the second step, we compute the elements of T>2 , that 
must be of the form A = Aidi + A2ii2 + X'j,d\ + A4(dic?2 + a%). 
Condition (i) yields A3 — A4 = 0, and together with (ii) we 
form the system 





1 -1 
1 -1 



A 4 



(11) 



J http:/ /www. neiu.edu/ zzeng/apatools.htm 



with columns indexed by d\,d\, d2,d\d2 + d|. We get two 
vectors in the kernel, the first yielding again di + d-2 and a 



second one for Ai = —1, A2 = 0, A3 = A4 = 1, so we deduce 
that — di + d\ + d\d 2 + d\ is a new element of T> 2 . 
In the third step we have 

A =Aidi + \ 2 d 2 + A 3 c!i + A 4 (did 2 + d\)+ 

As(di — d\) + A 6 (d2 + did 2 + d\d 2 — d\di), 

condition (i) leads to A3 — A4 + (A5 — Ae)di + (A5 — Ae)d2 = 0, 
and together with condition (ii) we arrive to 






-1 
-1 





" Ai " 




. A 6 _ 



(12) 



Since the kernel of this matrix gives elements that are al- 
ready in T> 2 , we derive that T> — T> 2 — T>z and the algorithm 
terminates. 

Note that for this example Macaulay's method ends with 
a matrix of size 12 x 10, instead of 4 x 6 in this approach. □ 

3.3 Computing a primal-dual pair 

In this section we provide a process that allows simulta- 
neous computation of a basis pair (D,B) of & and R/Q. 

Computing a basis of 3> degree by degree involves du- 
plicated computations. The successive spaces computed are 
Si C •■■ C 3>n = S>n+i- It is more efficient to compute only 
new elements AG ^ which are independent in 3lt/2>t-i at 
step t. 

Also, once dual basis is computed, one has to transform it 
into the form (|5). in order to identify a basis of R/Q as well. 
This transformation can be done a posteriori, by finding a 
sub-matrix of full rank and then performing Gauss- Jordan 
elimination over this sub-matrix, to reach matrix form ©. 

We introduce a condition (iii) extending Th. 13.21 that ad- 
dresses these two issues: It allows the computation of a total 
of n independent elements throughout execution, and re- 
turns a "triangular basis", e.g. a basis of R/Q is identified. 



Lemma 3.4. Let V t -i = (Ai, . 
whose coefficient matrix is 



,Ak) be a basis of S>t- 



Pi 
1 










»k 71 



1 



(13) 



yielding the standard basis B = (a? )i=i,...,fc. An element 
A 6 K[d] is not zero in @t/@t-i iff in addition to (i), (ii) 
of Th. \3.2\ we impose: 

(m) A[x Pi ] = 0, 1 < i < k ■ 

Proof. Let A £ K[d] be a non-zero functional satisfying 
(i), (ii) and (iii). Then A 6 @t and A[x 0i ] = for i = 
1, ...,k. If A £ St-u then A = £* =1 A t A t . Take for 
io the minimal i such that A^ 7^ 0. Then Afcc^'o] = A; , 
which is in contradiction with condition (iii). Thus, the non- 
zero solutions of (i), (ii) and (iii) correspond to the elements 
which are not zero in @t/@t-i- □ 

The above constraint is easy to realize; it is equivalent 
to Vz, d^ 1 ^ suppA 1 ^, which implies adding a row (linear 
constraint) for every i. In many cases this constraint is just 
Aifc = for some i, k, thus we rather remove the column 



corresponding to A« instead of adding a row. Either way, 
this lemma allows to shrink the kernel of the matrix and 
compute only new dual elements. 

Let us explore our running example, to demonstrate the 
essence of this improvement. 



Example 3.5. We re-run Ex. 13.31 using Lem. 13.41 
In the initialization step T>o — (1) is already in triangular 
form with respect to Bo = {1}. For the first step, we demand 
A[l] = 0, thus the matrix is the same as (|10p . yielding 
T>\ = (1, di + d 2 ). We extend B\ = {1,2:2}, so that T)\ is 
triangular with respect to B\. 

In the second step we remove from (|11|) the second colu- 
mn, hence 





" 1 -1 " 




Ai 






1 1 




A 3 


= 0, 




1 1 




A 4 




yielding a single solution — d\ H 


- d\ + did 2 + d\. We extend 



Bi by adding monomial Xi: B\ = {1, a; 2 , 21}. 

For the final ste p, we search an element with A[xi] 
two columns: 



A[x 2 ] = thus {12 



1- 
1-10 
10-10 
10 





" A3 - 




. A 6 . 



We find an empty kernel, thus we recover the triangular 
basis T> = T>2, which is then diagonalized to reach the form: 



Ai 
A2 
A3 



1 d 2 di d\ d\d 2 d\ 

1 

10 1 1 1 

1-1-1-1 



This diagonal basis is dual to the basis B = (1,0:2, xi) of 
the quotient ring and also provides a normal form algorithm 
(Lem. 12. 4p wrt B. In the final step we generated a 4 x 4 
matrix, size smaller compared to all previous methods. □ 

This technique for computing B can be applied similarly to 
other the matrix methods, e.g. Macaulay's dialytic method. 

If h(t) — h{t — 1) > 1, ie. there are more than one ele- 
ments in step t, then the choice of monomials to add to B is 
obtained by extracting a non-zero maximal minor from the 
coefficient matrix in (d a ). In practice, we will look first at 
the monomials of smallest degree. 

3.4 Approximate dual basis 

In our deflation method, we assume that the multiple 
point is known approximately and we use implicitly Taylor's 
expansion of the polynomials at this approximate point to 
deduce the dual basis, applying the algorithm of the previ- 
ous section. To handle safely the numerical problems which 
may occur, we utilize the following techniques: 

• At each step, the solutions of linear system i-iii) 
are computed via Singular Value Decomposition. Using a 
given threshold, we determine the numerical rank and an 
orthogonal basis of the solutions from the last singular values 
and the last columns of the right factor of the SVD. 

• For the computation of the monomials which define the 
equations (|3.4l iii) at the next step, we apply QR decom- 
position on the transpose of the basis to extract a non-zero 
maximal minor. The monomials indexing this minor are 
used to determine constraints ((5] i-iii). A similar numerical 
technique is employed in [21| . for Macaulay's method. 



4. DEFLATION OF A SINGULAR POINT 

We consider a system of equations / = (/i, . . . , f s ), which 
has a multiple point at x — £. Also, let B = (61, . . . , b^) be 
a basis of R/Q( and 2? = (Ai, . . . , A M ) its dual counterpart, 
with Ai = 1. 

We introduce a new set of equations starting from /, as 
follows: add for every fi the polynomial g k = f k + p k , p k = 
J2i=i £ i,A where e k = (e t>k , . . . , e Mife ) is a new vector of (i 
variables. 

Consider the system 



Vg(x,e) = [A^d^lg 



,A M (0.)to] 

where A 05 [<?*,] = Ai(9as)[fffc] is defined as in with £ re- 
placed by a;, ie. we differentiate g k but we do not evaluate 
at £. This is a system of (is equations, which we shall index 
T>g(x,e) = (51,1, . . - ,g^,s). We have 



g ik {x,e) = Af[f k +p k ] 



= A*[f k ]+A*\p k ]=A?[f k ]+p ik (x,e) 
A^[p k ] = e iik because V = (Ai, .., A M 



Notice that pi,fc(£, e) 
is dual to 6. 

As the first basis element of I? is 1 (the evaluation at the 
root), the first s equations are g(x,e) = 0. 

Note that this system is under-determined, since the num- 
ber of variables is /is + n and the number of equations is 
(is. We shall provide a systematic way to choose n variables 
and purge them (or better, set them equal to zero). 

This way we arrive to a square system T>g(x, i) (we use i 
for the remaining (is — n variables) of size (is x (is. We 
shall prove that this system vanishes on (£, 0) and that 

J D9 (C,o)^o. 

By linearity of the Jacobian matrix we have 

J-D B {x,e) = J-Df{x,e) + J-D P {x,e) 

= [Ji>f(x)\0 ] + [J% p (x, e)\J% p (x,e)], 

where J% p {x, e) (resp. J% p (x, e)) is the Jacobian matrix of 
Dp with respect to x (resp. e). 

Lemma 4.1. The Jacobian J% p (x,e) of the linear system 
~Dp = (pi,i, ■ ■ ■ ,Pn,s) withp iik {e k ) = Af [p k ](x, e k ) evaluated 
at (x,e) = (C, 0) is the identity matrix in dimension (is. 

Proof of Lemma 14.11 First note that the system is block 
separated, i.e. every p ik depends only on variables £; and not 
on all variables e = (ei, . . . ,£ n ). This shows that J p (x,e) 
is block diagonal, 



J% p {x,e) 



Ji 







Now we claim that all these blocks are all equal to the 
identity matrix. To see this, consider their entry d £kj \pi k ] 
for i,j = l,...,n, which is 



/-AtiW = Af[-^-p k ] = Ailbi] = 



, otherwise 



SmCe dJJjPk 



i,k) 



□ 



Lemma 4.2. The (is x n Jacobian matrix Jx>f{x) of the 
system TJf(x) = (/1, . . . , ff±n) is of full rank n at x = £. 

PROOF of Lemma 14.21 Suppose that the matrix is rank- 
deficient. Then there is a non-trivial vector in its kernel, 

J- Df (C)-v = 0. 



for i — 1, . . . . 



The entries of v are indexed by di. This implies that a non- 
zero differential A = vidi + ■ ■ ■ + v n d n of order one satisfies 
the following relations: (AAj)^[/j] = 0,i = l,...,(i,j = 
1, . . . , s. By the standard derivation rules, we have 

-nr ( AAi) = v k A t + A— A», 
dd k dd k 

fjt, , k = 1, . . . , n. Since 3) is stable by deriva- 
3. We deduce that the vector space spanned 
by {2>, AS>) is stable by derivation and vanishes on / at £. 
By Proposition E2 we deduce that C 3>. This is a 
contradiction, since A is of degree 1 and the the elements in 
3) are of degree < N. □ 

The columns of Jx> g (x, e) are indexed by the variables (a;, e), 
while the rows are indexed by the polynomials g ik . We con- 
struct the following systems: 

(a) Let £>/ 7 be a subsystem of T>f such that the corre- 
sponding n rows of J-Df(C) are linearly independent 
(Lem. g?2]). We denote by I = {(h,ki), . . . , (i n ,k n )} 
their indices. 

(b) Let T)g(x, e) be the square system formed by removing 
the variables £i ltkl , £i n ,k n iroxa T>g(x,e). There- 
fore the Jacobian Jx>g(x, e) derives from J-p g (x, e), af- 
ter purging the columns indexed by Si 1:kl , . . . ,£i„, kn , 



and it's (ij,kj) row becomes [V(Af. §i 







Theorem 4.3 (Deflation Theorem 1). Let f(x) be 
a n—variate polynomial system with an (i—fold isolated zero 
at x — £. Then the n x n system Vf I (x) = 0, defined in 

(a) , has a simple root at x — £. 

Proof. By construction, £ is a solution of Vf I (x) — 0. 
Moreover, the indices I are chosen such that det Jx>fi(C) 7^ 
0. This shows that £ is a simple (thus isolated) root of the 
system Vf'ix) = 0. □ 

Example 4.4. In our running example, we expand the 
rectangular Jacobian matrix of 6 polynomials in (xi,X2). 
Choosing the rows corresponding to /1 and (di — d\ — d\di — 
we nn d a non-singular minor, hence the resulting 
system (/i,2xi) has a regular root at C, = (0,0). □ 

The deflated system Df 1 (x) = is a square system in n 
variables. Contrarily to the deflation approach in [7] [3] , we 
do not need to introduce new variables here and one step 
of deflation is sufficient. In the following theorem, we do 
introduce new variables to express the condition that the 
perturbed system has a given multiplicity structure. 

Theorem 4.5 (Deflation Theorem 2). Let f(x) be 
a n~variate polynomial system with an p—fold isolated zero 
at x — £. The square system T>g(x,e) = 0, as defined in 

(b) , has a regular isolated root at (x, e) — (£, 0). 

Proof. By definition of V, we have 

Vg(C,0) = (A< [/],..., A«[/])=0. 

Moreover, by construction of T>g we get, up to a row per- 
mutation, the determinant: 



±det J-Dg(CO) 



det 



Ji 



det Ji / 0, 



where J\ = Jx>f(C)- This shows that (£, 0) is regular 
and thus isolated point of the algebraic variety defined by 
Vg(x,e) = 0. □ 



Nevertheless, this deflation does differ from the deflation 
strategy in [7] [3]. There, new variables are added that cor- 
respond to coefficients of differential elements, thus introdu- 
cing a perturbation in the approximate dual basis, in case 
of approximate input. Hence the output concerns a deflated 
root of the given approximate system. In our method, we 
perturb the equations, keeping an approximate structure of 
a multiple point. Consequently, the certification of a root 
concerns a nearby system, within controlled error bounds, 
as it shall be described in Sect.JS] 

We mention that it would also be possible to use the equa- 
tions ^ i-iii) to construct a deflated system on the differ- 
entials and to perturb the approximate dual structure. 

5. VERIFYING APPROXIMATE SINGULAR 
POINTS 

In real-life applications it is common to work with ap- 
proximate inputs. Also, there is the need to (numerically) 
decide if an (approximate) system possesses a single (real) 
root in a given domain, notably for use in subdivision-based 
algorithms, eg. [Hfl [10] . 

In the regular case, Smale's a— theory, extending New- 
ton's method, can be used to answer this problem. Another 
option is Rump's Theorem, also based on Newton theory. 
In our implementation we choose this latter approach, since 
it is suitable for inexact data and suits best with the pertur- 
bation which is applied. Our perturbation coincides to the 
numerical scheme of [16] in the univariate case. 

The certification test is based on the verification method 
of Rump |16l Th. 2.1], which we rewrite in our setting: 

Theorem 5.1 ([H] Rump's Theorem). Let f £ R n 
be a polynomial system and £* £ R" a real point. Given 
an interval domain Z £ IR™ containing £* £ R n , and an 
interval matrix M £ IR nxn whose i—th column Mi satisfies 
Vfi(Z) Mi for i — 1 . . . , n, then the following holds: 

If the inclusion 

v(f, z, o = -jficr'fin + j f (cr 1 M)z ci 

(14) 

is true, then there is a unique £ £ Z with / (C) = and the 
Jacobian matrix Jf (£) £ M is non-singular. 

This theorem is applied on the perturbed system. If the 
test succeeds, we also get a domain for £— variables that re- 
flects the distance of the approximate system from a precise 
system with the computed local structure. 

Example 5.2. We start with an approximate system: f\ = 
1.071xi-1.069a;2 + 1.018a; 1 , f 2 = 1.024a;i-1.016a;2 + 1.058a;| 
and the domain: Z = [-.01, .03] x [-.03, .01]. The Jacobian 
at x = (0, 0) evaluates to .00652, hence it is close to singular. 

We set the tolerance equal to .04, i.e. the size of the 
domain, and we consider the center of the box as our appro- 
ximate point, £* = (.01,— .01). 

First we compute approximate multiplicity structure at 
C*, V = (1, d 2 + 1.00016di + .99542^2 + 1.03023d?, di - 
1.00492d| - 1.00016did 2 - 1.03514d?) as well as (l,x 2 ,x 1 ), 
a standard basis for the quotient. The latter indicates to 
perturb up to linear terms. 

Now apply the second deflation theorem 14.51 to get the 
6x6 system g = (1.018 x\ + 1.071 an + (e 12 - 1.069) x 2 , e 12 - 
.02023, .01723 + 2.036x1, 1.058 x\ + (1.024 + £23)2:1 + (£22 - 



1.016)2:2 +£21, .04217 + 2.116 2:2 +£22, £23- .03921), which 
has a regular root for £ £ Z and parameters (£12, £21 , £22, £23). 
Indeed, applying Theorem Owith Z' = Zx [-.04, .04] 4 and 

O 

(C*,0, ..,0) we get an inclusion V(g,Z',C) QZ' . □ 

6. GEOMETRY AROUND A SINGULARITY 

As a final step in analyzing isolated singularities, we show 
how the local basis can be used to compute the topological 
degree around the singular point. If the latter is a self- 
intersection point of a real algebraic curve, one can deduce 
the number of curve branches that pass through it. 

Topological degree computation. Let f(x) be a square 
n— variate system with an fi— fold isolated zero at x = £. To 
a linear form A £ R[9<;], we associate the quadratic form 

Qa ■ R/Q x E/Q -> R , (bi,bj) h-> A(bibj) (15) 

for R/Q = (6i,...,6 M ). The signature of this (symmetric 
and bi-linear) form is the sum of signs of the diagonal entries 
of any diagonal matrix representation of it. 

Proposition 6.1 (d Th. 1.2]). IfQ$, $ £ V is any 
bi-linear symmetric form such that $^[det Jf{x)\ > 0, then 

tdeg 6 (f) = sgn(Qs). (16) 

This signature is independent of the bi-linear form used. 

We can use this result to compute the topological degree 
at x = £ using the dual structure at £. Since a basis T> is 
available we set $ = ± Ai , for some basis element that is not 
zero on det Jf(x). Indeed, such an element can be retrieved 
among the basis elements, since det Jf ^ (/), see [5] Ch. 0]. 

In practice it suffices to generate a random element of 
compute it's matrix representation [&(bibj)]ij , and then 
extract the signature of Q$ . 

Branches around a singularity. In the context of com- 
puting with real algebraic curves, the identification of self- 
intersection points is only the first step of determining the 
local topology. As a second step, one needs to calculate the 
number of branches attached to the singular point £, here- 
after denoted Br(/,£). This information is encoded in the 
topological degree. 

An implicit curve in n— space is given by all points sati- 
sfying f(x) = 0, / = (/1, . . . , /n-i). Consider p(x) = (xi - 
Ci) 2 + ■ • ■ + (x n - <n) 2 , and 5 (a;) = det J (M (as). Then ([E] 
and references therein): 

Br(/,C) = 2tdeg f (/, 5 ). (17) 

This implies an algorithm for Br(/,£). First compute the 
primal-dual structure of (/,<?) at £ and then use Prop. 16.11 
to get tdeg c (/,#). 

Example 6.2. Consider the implicit curve f(x, y) = 0, in 
xy~ plane, with f(x, y) = x 4 + 2x 2 y 2 + y 4 + 3x 2 y — y 3 , that 
looks like this We search for the number of branches 

touching C = (0,0). 

We compute g(x,y) = J^f x 2 +y 2~ j = l&xy 2 — 6x 3 , and then 
the multiplicity structure of (/, g) at £, and we arrive to the 
standard basis B = (l,y,x,y 2 ,xy,x 2 ,y 3 ,xy 2 ,x 2 y). Among 
the 9 elements of the dual basis, we find $ = 0% + |d 4 + 
\d 2 d 2 + |d 4 , having value $°[det J(f }9 ) (x)] = 54 > on the 
Jacobian determinant. 



Using B and ([6]), we get the 9x9 matrix representation of 
Q<s> (fl5)) with ij— th entry $[NF(cc' 3i a;' 3 J ')], and we compute 
tdeg c (/, g) = sgn = 3, thus Br(/, (0, 0)) = 6. □ 

7. EXPERIMENTATION 

Our method is developed in Maple. It uses Mourrain's 
Integration technique to compute (approximate) dual basis 
and derive the augmented system of Th. 14.51 Then Rump's 
method is used to verify the root. Macaulay's method is also 
implemented for testing purposes. 

Example 7.1. Consider the system [S] of 3 equations in 2 
variables fi=xl + xix%, h = X\x% + x%, fa = xlx^+xixl, 
and the singular point (0,0). 

Suppose that the point is given. Using 13.21 and 13.41 we 
derive the primal-dual pair T> = (1, di, (I2, d\, did%, d\, d,2 + 
dl + d\ d2 — di d\ ) , where d| is underlined to show that it cor- 
responds to x\ in the primal standard basis B. The biggest 
matrix used, in depth 4, was of size 9x8, while Macaulay's 
method terminates with a matrix of size 30 x 15. 

To deflate the root, we construct the augmented system 
T>f of 21 equations. The 21 x 2 Jacobian matrix Jf;^ is of 
rank 2 and a full-rank minor consists of the rows 4 and 5. 
Consequently find the system di<h[fi]) = (3xi,2x2) 

which deflates (0, 0). Note that even though both equations 
of the deflated system derive from ft, the functionals used 
on fi are computed using all initial equations. □ 

Example 7.2. Let, as in [B] |S], /1 = 2a: 1 + 2x\ + 2x 2 + 
2x2 + oos — 1, /2 = (xi + x 2 - x z - l) 3 - x\, and f 3 = 
2a;? + 2x1 + Wx 3 + 5a;§ + 5) 3 - 1000a;?. 

Point (0,0,-1) occurs with multiplicity equal to 18, in 
depth 7. The final matrix size with our method is 206 x 45, 
while Macaulay's method ends with a 360 x 165 matrix. 

If the objective is to deflate as efficiently as possible, then 
one can go step by step: First compute a basis for 5?i and 
stop the process. We get the evaluation 1 and 2 first order 
functionals, which we apply to /1. We arrive to (l[/i], (efo — 
di)[.fi], (di+rf 3 )[/i]) = {fi,-Axi+4x2,2 + 4 Xl +2x3) and 
we check that the Jacobian determinant is 64, thus we have 
a deflated system only with a partial local structure. □ 

Table [7] shows computations on the benchmark set of [3]- 
Multiplicity, matrix sizes at termination step are reported. 



Sys. 




Integration 


Macaulay 


cmbsl 


11 


33 x 23 


105 x 56 


cmbs2 


8 


21 x 17 


60 x 35 


mthl91 


4 


10 x 9 


30 x 20 


decker2 


4 


5x5 


20 x 15 


Ojika2 


2 


6x5 


12 x 10 


Ojika3 


4 


24 x 9 


60 x 35 


KSS 


16 


569 x 69 


630 x 252 


Caprasse 


4 


34 x 13 


60 x 35 


Cyclic 9 


4 


369 x 33 


495 x 33 


DZ1 


131 


1450 x 524 


4004 x 1365 


DZ2 


7 


73 x 33 


360 x 165 


DZ3 


5 


14 x 8 


30 x 21 



Table 1: Benchmark systems from [3]. 
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APPENDIX 

We attach additional examples that did not fit page limits. 
Verification Example. 

Let /i = (x\x2 — X\x%,f2 = xi — 22). The verification 
method of [TB] applies a linear perturbation on this system, 
but fails to certify the root x — (0, 0). 

We consider an approximate point £* = (.01, .002) and 
we compute the approximate multiplicity structure: 

V = (Ai,...,A 4 ) = (1.0, l.Ocfe, l.Orfi +l.Orfj, LOrfida +l.Odjj) 

The augmented system g(x) — (A 3 (/i)) = {fx, 2.02:12:2 — 
1.02^-1.02:1, 2.O21-2.O22, I.O21-I.O2I, / 2 , -2.O22, 0.,0.) 
has a Jacobian matrix: 



.00 
.00 



.016 -.99 
-.02 .016 



2.0 
-2.0 



1.0 

-.004 




-2.0 



with a non-zero minor at the third and forth row. Using 
this information, we apply the following perturbation to the 
original system: 

fx = x\x2 - x x x\ + en + £l 2 2 2 

fs = 2l - 2a + £21 + £2222 + £2321 + £242l2 2 

Thus g(xi, 22, £11, £12, £21, £22, £23, £24), computed as before, 
is a square system with additional equations 



h 
h 
U 
h 
h 

fs 



1.02? 



2.O2122 + l.Oeia 



I.O21 



2.O2122 — 1.02 2 
2.O21 - 2.O22 

-2.022 + 1.0£ 22 + 1.02i£24 

l-0£23 + 1.022E24 

1.0£24 



Now take the box Zi = [-.03, .05] x [-.04, .04] x [-.01, .01] 6 . 
We apply Th. 15. li on g, ie. we compute V(g, Z\ , £* ) . For the 
variable £21 the interval is [—.015, .15] % (—.01, .01), there- 
fore we don't get an answer. 

We shrink a little Z x down to Z 2 = [-.03, .05] x [-.02, .02] x 
[-.01, .01] 6 and we apply again Th. 15.11 which results in 

[-.004, .004] 
[-.004, .004] 
[-.001, .001] 
[-.007, .007] 
[-.006, .006] 
[-.009, .009] 
-.00045, .00035] 
[.0, .0] 

thus we certify the multiple root of the original system inside 
Z 2 . □ 



V(g,Z 2 ,C) = 



(=Z2, 



Example 16.21 (Cont'd). 

Consider the implicit curve f(x, y) = 0, in xy— plane, with 
f(x, y) = 2 4 + 2x 2 y 2 +y i + 3x 2 y - y :i , 

that looks like this We search for the number of branches 
touching < = (0,0). 

This point is of multiplicity 4, as the dual basis we get for 

is (1, d x , d y , d 2 + dy), which provides no information for the 
number of branches. 
We compute 

> 



g(x,y) = J (fiX 2 +y 2 ) = 18xy 

and then the multiplicity structure of (/, 
arrive to 



6z 3 , 



at £, and we 



D — (1, dy, d x , d yl d x dy , d x , d y + d y -(- d x d y + d x , 

— OO o 

, ,2 1 , ,3 ,2 , 9 ,4 _ 3 2 j2 _ 9 4 . 

U> X CLy + O flj,, Uj x CLy dy U, x U,y CL X J, 

00 o 

and the standard basis 

B = (l,y,x,y 2 ,xy,x 2 ,y 3 ,xy 2 ,x 2 y). 
Among the 9 elements of the dual basis, we find 

* = d y + ^d y + ^d 2 x d 2 y + ^d^, 

having value $°[det J(j s )(a;)] = 54 > on the Jacobian 
determinant. 

Using B and ([6]), we compute the matrix NF(x l3i x l3 i ) of 
multiplication in R/ (/, g), given at the end of the page. Now 
a representation of Qa> (|15[) can be computed, by applying 
$° on the multiplication table to get: 
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With a QR iteration, we find 6 positive and 3 negative 
eigenvalues of this representation, hence we compute 

tdeg c (/, g) = sgn Q e = 6 - 3 = 3, 

i.e. there are 6 branches of the curve around (0, 0) . □ 



1 


y 


X 


V 2 


xy 
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y 3 
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x 2 y 
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xy 
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x 2 y 







\y 3 - l x2 y 
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xy 


.v 2 
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xy 


x 2 y 


3xy 2 





l/8y 3 -3/8x 2 y 





y 2 


y 3 


xy 2 







- h 2 y 











xy 
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xy 
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x £ y 





W - h 2 y 














X 2 




3xy 2 


\y 3 - 3/8x 2 j/ 





3/8^3 _ 9 x 2 y 











y 3 


3/8y 3 - Ix^y 























xy z 





§2/ 3 - lx 2 y 
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x y 


l/8j,3 - §x2<y 
























Multiplication table for R/(f,g) (Example |6T2|) . 



